#include <zMat.hpp>
#include <UnitTest++.h>

#ifdef ZZZ_LIB_LAPACK
SUITE(zMatrixMathTest)
{
  TEST(LU)
  {
    zzz::zMatrix<double> x(3,3);
    x(0,0)=1;
    x(0,1)=2;
    x(0,2)=3;
    x(1,0)=3;
    x(1,1)=2;
    x(1,2)=1;
    x(2,0)=3;
    x(2,1)=1;
    x(2,2)=2;
    zzz::zMatrix<long> IPIV;
    CHECK(LUFactorization(x,IPIV));
    CHECK_CLOSE(3.0000,x(0,0),0.0001);
    CHECK_CLOSE(2.0000,x(0,1),0.0001);
    CHECK_CLOSE(1.0000,x(0,2),0.0001);
    CHECK_CLOSE(0.3333,x(1,0),0.0001);
    CHECK_CLOSE(1.3333,x(1,1),0.0001);
    CHECK_CLOSE(2.6667,x(1,2),0.0001);
    CHECK_CLOSE(1.0000,x(2,0),0.0001);
    CHECK_CLOSE(-0.7500,x(2,1),0.0001);
    CHECK_CLOSE(3.0000,x(2,2),0.0001);
  }

  TEST(INVERT)
  {
    zzz::zMatrix<double> x(3,3);
    x(0,0)=1;
    x(0,1)=2;
    x(0,2)=3;
    x(1,0)=3;
    x(1,1)=2;
    x(1,2)=1;
    x(2,0)=3;
    x(2,1)=1;
    x(2,2)=2;
    zzz::zMatrix<double> ix(x);
    Invert(ix);
    CHECK_CLOSE(-0.2500,ix(0,0),0.0001);
    CHECK_CLOSE(0.0833,ix(0,1),0.0001);
    CHECK_CLOSE(0.3333,ix(0,2),0.0001);
    CHECK_CLOSE(0.2500,ix(1,0),0.0001);
    CHECK_CLOSE(0.5833,ix(1,1),0.0001);
    CHECK_CLOSE(-0.6667,ix(1,2),0.0001);
    CHECK_CLOSE(0.2500,ix(2,0),0.0001);
    CHECK_CLOSE(-0.4167,ix(2,1),0.0001);
    CHECK_CLOSE(0.3333,ix(2,2),0.0001);
    zzz::zMatrix<double> a(x*ix);
    CHECK_CLOSE(1,a(0,0),0.0001);
    CHECK_CLOSE(0,a(0,1),0.0001);
    CHECK_CLOSE(0,a(0,2),0.0001);
    CHECK_CLOSE(0,a(1,0),0.0001);
    CHECK_CLOSE(1,a(1,1),0.0001);
    CHECK_CLOSE(0,a(1,2),0.0001);
    CHECK_CLOSE(0,a(2,0),0.0001);
    CHECK_CLOSE(0,a(2,1),0.0001);
    CHECK_CLOSE(1,a(2,2),0.0001);
  }
  TEST(SVD)
  {
    zzz::zMatrix<double> x(3,3);
    x(0,0)=1;
    x(0,1)=2;
    x(0,2)=3;
    x(1,0)=3;
    x(1,1)=2;
    x(1,2)=1;
    x(2,0)=3;
    x(2,1)=1;
    x(2,2)=2;
    zzz::zMatrix<double> U,S,VT;
    zzz::zMatrix<double> x1(x);
    CHECK(SVD(x1,U,S,VT));
    x1=U*zzz::Diag(S)*VT;
    for (zzz::zuint r=0; r<x.Size(0); r++) for (zzz::zuint c=0; c<x.Size(1); c++)
      CHECK_CLOSE(x(r,c),x1(r,c),zzz::EPSILON);
  }

  TEST(DET)
  {
    zzz::zMatrix<double> x(3,3);
    x(0,0)=1;
    x(0,1)=2;
    x(0,2)=3;
    x(1,0)=3;
    x(1,1)=2;
    x(1,2)=1;
    x(2,0)=3;
    x(2,1)=1;
    x(2,2)=2;
    double d=Det(x);
    CHECK_CLOSE(-12,d,zzz::EPSILON);
  }

  TEST(EIGEN)
  {
    zzz::zMatrix<double> x(3,3);
    x(0,0)=1;
    x(0,1)=2;
    x(0,2)=3;
    x(1,0)=2;
    x(1,1)=4;
    x(1,2)=5;
    x(2,0)=3;
    x(2,1)=5;
    x(2,2)=6;
    zzz::zMatrix<double> eigen;
    CHECK(EigenSym(x,eigen));

    CHECK_CLOSE(0.7370,x(0,0),0.0001);
    CHECK_CLOSE(0.5910,x(0,1),0.0001);
    CHECK_CLOSE(0.3280,x(0,2),0.0001);
    CHECK_CLOSE(0.3280,x(1,0),0.0001);
    CHECK_CLOSE(-0.7370,x(1,1),0.0001);
    CHECK_CLOSE(0.5910,x(1,2),0.0001);
    CHECK_CLOSE(-0.5910,x(2,0),0.0001);
    CHECK_CLOSE(0.3280,x(2,1),0.0001);
    CHECK_CLOSE(0.7370,x(2,2),0.0001);
    CHECK_CLOSE(-0.5157,eigen[0],0.0001);
    CHECK_CLOSE(0.1709,eigen[1],0.0001);
    CHECK_CLOSE(11.3448,eigen[2],0.0001);
  }
}
#endif // ZZZ_LIB_LAPACK